1 Abstract

Most data has a spatial dimension to it. Knowing ‘where’ the data is coming from is often as crucial as knowing the ‘what’, ‘when’ and ‘who’ dimensions of a given dataset. Therefore it should be no surprise that R has a rich suite of packages for constructing maps and analyzing spatial data. R’s capability has grown so much over the years that it’s functionality rivals many dedicated geographic information systems (GIS). During this Meetup the basics for managing and mapping spatial data with be introduced, using the following packages: sf, ggplot2, tmap, mapview and leaflet.

2 Example datasets

How you plot your data in R depends on what format it’s in. R has several different formats for managing spatial data (e.g. sp vs sf), and different formats are only compatible with certain plotting systems. Thankfully converting between different spatial formats is not difficult.

2.1 maps R package

library(maps)
library(maptools)

map(database = "state")

map(database = "state", region = "indiana")

st <- map("state", fill = TRUE, plot = FALSE)
st_sp <- map2SpatialPolygons(st, IDs = st$names)
st_sp$state <- st$names
proj4string(st_sp) <- CRS("+init=epsg:4326")

# data frame
str(st)
## List of 4
##  $ x    : num [1:15599] -87.5 -87.5 -87.5 -87.5 -87.6 ...
##  $ y    : num [1:15599] 30.4 30.4 30.4 30.3 30.3 ...
##  $ range: num [1:4] -124.7 -67 25.1 49.4
##  $ names: chr [1:63] "alabama" "arizona" "arkansas" "california" ...
##  - attr(*, "class")= chr "map"
# sp object
str(st_sp, 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  63 obs. of  1 variable:
##   ..@ polygons   :List of 63
##   ..@ plotOrder  : int [1:63] 50 28 4 33 30 2 5 63 44 25 ...
##   ..@ bbox       : num [1:2, 1:2] -124.7 25.1 -67 49.4
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
plot(st_sp)

2.2 County Boundaries

library(sf)
library(USAboundaries)

cnty <- us_counties() 
cnty <- subset(cnty, !state_name %in% c("Alaska", "Hawaii"))
str(cnty)
## Classes 'sf' and 'data.frame':   3186 obs. of  13 variables:
##  $ statefp          : chr  "39" "46" "55" "72" ...
##  $ countyfp         : chr  "131" "003" "035" "145" ...
##  $ countyns         : chr  "01074078" "01266983" "01581077" "01804553" ...
##  $ affgeoid         : chr  "0500000US39131" "0500000US46003" "0500000US55035" "0500000US72145" ...
##  $ geoid            : chr  "39131" "46003" "55035" "72145" ...
##  $ name             : chr  "Pike" "Aurora" "Eau Claire" "Vega Baja" ...
##  $ lsad             : chr  "06" "06" "06" "13" ...
##  $ aland            : num  1.14e+09 1.83e+09 1.65e+09 1.19e+08 1.72e+09 ...
##  $ awater           : num  9567612 11201379 18848512 57805868 1496797 ...
##  $ state_name       : chr  "Ohio" "South Dakota" "Wisconsin" "Puerto Rico" ...
##  $ state_abbr       : chr  "OH" "SD" "WI" "PR" ...
##  $ jurisdiction_type: chr  "state" "state" "state" "territory" ...
##  $ geometry         :sfc_MULTIPOLYGON of length 3186; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:9, 1:2] -83.4 -82.8 -82.8 -83 -83.2 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  NA NA NA NA ...
vars <- c("statefp", "countyfp")
plot(cnty[vars])

2.3 USDA-NASS Corn Yield Data

https://github.com/potterzot/rnassqs

# devtools::install_github('potterzot/rnassqs')
# library(rnassqs)

# source("C:/Users/steph/Nextcloud/code/api_keys.R")

data(state)
st <- state.abb

# corn_us <- lapply(st, function(x) {
#   cat("getting", x, as.character(Sys.time()), "\n")
#   tryCatch({
#     corn = nassqs_yield(
#       list("commodity_desc"="CORN", 
#            "agg_level_desc"="COUNTY", 
#            "state_alpha"=x
#            ),
#       key = nass_key
#       )}, 
#     error = function(err) {
#       print(paste("Error occured:  ",err))
#       return(NULL)
#       }
#     )
#   })
# corn_us <- do.call("rbind", corn_us)
# 
# save(corn_us, file = "C:/workspace2/corn_us.RData")
# write.csv(corn_us, file = "nass_corn_us.csv", row.names = FALSE)

load(file = "C:/Users/Stephen.Roecker/Nextcloud/data/corn_us.RData")
corn_yield <- subset(corn_us, short_desc == "CORN, GRAIN - YIELD, MEASURED IN BU / ACRE")
corn_yield <- within(corn_yield, {
  Value      = as.numeric(Value)
  year       = as.numeric(year)
  state_name = NULL
  state      = state_alpha 
  })


cnty_corn <- merge(cnty, corn_yield, 
                   by.x = c("state_abbr", "countyfp"), 
                   by.y = c("state_alpha", "county_code"), 
                   all.x = TRUE
                   )

corn_states <- c("IL", "IA", "IN", "MI", "MN", "MO", "NE", "OH", "SD", "ND", "WI")

library(dplyr)
library(ggplot2)

group_by(corn_yield, state_alpha, year) %>%
  summarize(
    yield_low    = min(Value, na.rm = TRUE),
    yield_median = median(Value, na.rm = TRUE),
    yield_max    = max(Value, na.rm = TRUE)
    ) %>%
  filter(state_alpha %in% corn_states) %>%
  mutate(state = state_alpha,
         source = "NASS"
         ) %>%
  ggplot() +
  geom_line(aes(x = year, y = yield_median, col = source)) +
  geom_ribbon(aes(x = year, ymin = yield_low, ymax = yield_max), alpha = 0.25) +
  facet_wrap(~ state) +
  ylab("yield per county (bu/acre)") +
  ggtitle("USDA-NASS Corn Yields")

  # geom_point(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, col = "NASIS"), size = 1) +
  # geom_ribbon(data = yld_sum2[yld_sum2$state %in% corn_states, ], aes(x = year, ymin = yield_low2, ymax = yield_max2, col = "NASIS"), alpha = 0.25) +
  # geom_pointrange(data = yld_sum[yld_sum$state %in% corn_states, ], aes(x = 2018, y = yield_med, ymin = yield_low2, ymax = yield_max2, col = source))

2.4 Indiana General Soil Map

IN <- read_sf(dsn = "I:/geodata/soils/soils_GSMCLIP_mbr_2599033_03/wss_gsmsoil_IN_[2006-07-06]/spatial/gsmsoilmu_a_in.shp", layer = "gsmsoilmu_a_in")

# simplify polygons
IN <- rmapshaper::ms_simplify(IN)

2.5 Kellogg Soil Survey Laboratory Data

library(aqp)
library(soilDB)

# download lab locations for the Miami soil series
miami <- fetchKSSL("Miami")
idx <- complete.cases(miami$x, miami$y)
miami_sp <- miami[idx, ] # remove sites with missing locations

coordinates(miami_sp) <- ~ x + y
proj4string(miami_sp) <- CRS("+init=epsg:4326")

miami_sp <- SpatialPointsDataFrame(
  coords = coordinates(miami_sp@sp),
  data   = site(miami_sp),
  proj4string = CRS(proj4string(miami_sp))
  )
miami_sf <- st_as_sf(miami_sp)
miami_sf <- within(miami_sf, {
  lon = st_coordinates(miami_sf)[, 1]
  lat = st_coordinates(miami_sf)[, 2]
  })

3 Example maps with several R packages

3.1 ggplot2

library(ggplot2)

# loads data from maps package
st <- map_data("state")

# tidy example data
st_tidy <- broom::tidy(st_sp, region = "state") 
IN_tidy <- broom::tidy(as(IN, "Spatial"), region = "MUSYM")

dim(st_sp)
## [1] 63  1
dim(st)
## [1] 15537     6
# Lines
ggplot() +
  geom_point(data = miami_sf, aes(x = lon, y = lat)) +
  geom_path(data = st, aes(x = long, y = lat, group = group)) +
  xlim(range(miami_sf$lon)) +
  ylim(range(miami_sf$lat)) +
  ggtitle("Location of Miami Lab Pedons")

# Polygons
ggplot() +
  geom_polygon(data = st_tidy, aes(x = long, y = lat, group = group, fill = id)) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  # remove legend
  guides(fill = FALSE)

# sf package

# Polygons
ggplot() +
  geom_sf(data = cnty, aes(fill = statefp, lty = NA)) +
  geom_sf(data = miami_sf) +
  coord_sf(crs = "+init=epsg:5070") +
  guides(fill = FALSE)

# Facets

test <- cnty_corn[cnty_corn$year %in% 2012:2017, ]

ggplot() +
  geom_sf(data = test, aes(fill = Value, lty = NA)) +
  scale_fill_viridis_c(na.value = "transparent") +
  facet_wrap(~ year) +
  geom_path(data = st, aes(x = long, y = lat, group = group)) +
  ggtitle(corn_yield$short_desc[1])

3.2 ggmap

library(ggmap)

# get ggmap
bb <- sf::st_bbox(IN)
bb <- make_bbox(lon = bb[c(3, 1)], lat = bb[c(2, 4)])
gmap <- get_map(bb, maptype = "terrain", source = "osm")


# Lines
ggmap(gmap) +
  geom_path(data = IN_tidy, aes(x = long, y = lat, group = group))

  # geom_sf(data = IN, fill = NA, inherit.aes = FALSE) +
  # guides(fill = FALSE)

# geom_sf() doesn't work with ggmp, their is a systematic shift https://github.com/r-spatial/sf/issues/336

3.3 tmap

library(tmap)

tm_shape(IN) + tm_polygons("MUSYM", border.col = NULL) +
  tm_shape(cnty) + tm_borders() +
  tm_shape(miami_sf) + tm_dots() +
  tm_legend(legend.outside = TRUE)

# interactive web map 
tmap_mode("view")

tm_basemap("OpenStreetMap") +
  tm_shape(IN) + tm_borders()

3.4 mapview

library(mapview)

cols <- RColorBrewer::brewer.pal(50, "Paired")

test <- mapview(IN, zcol = "MUSYM", lwd = 0, col.regions = cols) +
  mapview(cnty, type = "l")

test
# export to html
mapshot(test, url = "C:/workspace2/test.html", selfcontained = FALSE)

3.5 leaflet

library(leaflet)

test <- leaflet() %>%
  addProviderTiles("Esri.WorldImagery", group = "Imagery") %>%
  addPolygons(data = IN, fill = FALSE, color = "black", weight = 2)

test
# export to html
htmlwidgets::saveWidget(test, file = "C:/workspace2/test.html", selfcontained = FALSE)

4 Additional Reading

Healy, K., 2018. Data Visualization: a practical introduction. Princeton University Press. http://socviz.co/

Gimond, M., 2019. Intro to GIS and Spatial Analysis. https://mgimond.github.io/Spatial/

Hijmans, R.J., 2019. Spatial Data Science with R. https://rspatial.org/

Lovelace, R., J. Nowosad, and J. Muenchow, 2019. Geocomputation with R. CRC Press. https://bookdown.org/robinlovelace/geocompr/

Pebesma, E., and R. Bivand, 2019. Spatial Data Science. https://keen-swartz-3146c4.netlify.com/